# Function to prepare Data and run stan to sample from posteror 
######################

est_stan <- function(data,iteration=200, chains=1, country, year,
                     pos=paste("pos.",c(1:3,5:10,13,16),sep=""),
                     model="model_graded_response.stan",
                     ...){
  
  # Install required Packages
  require(reshape2)
  require(rstan)  
  rstan_options(auto_write = TRUE) # Multipel Cores
  options(mc.cores = parallel::detectCores())
  
  # Prepare Data
  sel  <- c("ctry_name","pid","eid",paste("pos.",c(1:3,5:10,13,16),sep=""))
  data <- data[data$ctry_name==country & data$year==year,]
  cd   <- melt(data[,sel], id.vars = c("ctry_name","pid","eid"))
  cd   <- na.omit(cd)
  
  # recode items 1:K
  cd$variable <- as.character( cd$variable)
  ui <- unique(cd$variable)
  for(i in 1:length(ui)) cd$variable[cd$variable==ui[i]] <- i
  
  # Stan Data
  uc <- function(x) length(unique(as.factor(x)))
  standata <- list(T = uc(cd$value), J = uc(cd$pid), K = uc(cd$variable),N = nrow(cd), 
                   I = uc(cd$eid),
                   jj = as.numeric(as.factor(cd$pid)),
                   kk = as.numeric(cd$variable),
                   ii = as.numeric(as.factor(cd$eid)),
                   y = cd$value +1 )
  
  # Initial Value Function; lstoetze: nu is fixed otherwise nummerical errors
  init_fun <- function(T = uc(cd$value), J = uc(cd$pid), K = uc(cd$variable),I = uc(cd$eid)) {
    list(
      tau = sort(runif(T-1,-2,2))
      ,bk_free  = runif(K,1,1) 
      ,nu_free = runif(J,1,1)
    )
  }
  
  # Call stan
  fit <- stan(file = model, data = standata, init = init_fun,
              iter = iteration, chains = chains, ...)
  
  # Return Fit
  return(fit)
}

